DADES D’EIVISSA I FORMENTERA
Maria Àngels Jaume
2023-03-10
1. LECTURA DE DADES
En aquest document farem el mateix estudi que amb les dades de Mallorca però utilitzant les d’Eivissa i Formentera. Els càlculs i procediments són anàlegs al de l’illa de Mallorca, per això obviarem alguns detalls.
En primer lloc, grafiquem en forma de sèrie temporal aquestes dades:
turisme <- read_excel("IBESTAT.xls")
turisme$Data <- gsub("M","-",turisme$Data)
gastos.ts<-ts(turisme[-1], start=c(2015,10), frequency = 12)
ef <- data.frame(x=1:86, y=turisme$`Eiv-Form`) #dades d'Eivissa i Formentera
serie_ef <- ts(ef$y,start=c(2015,10),frequency = 12)
plot.ts(serie_ef, main="Eivissa i Formentera", xlab="Any", ylab="Despeses mensuals en €")
Les dades van des d’octubre de 2015 a novembre de 2022 (llavors tenim 7 cicles complets). Podem observar que l’efecte de la COVID és veu clarament en 2020.
Per veure millor l’estacionalitat de la sèrie visualitzem cada període mensualment:
seasonplot(serie_ef, col=c("brown", "blue","red", "orange", "pink", "purple", "yellow","green"),year.labels=TRUE, main="Estacionalitat d'Eivissa i Formentera", xlab="Mes", ylab=" Despeses en €")
legend("bottomright",
legend = c(2015, 2016, 2017, 2018, 2019, 2020, 2021, 2022),
fill = c("brown", "blue","red", "orange", "pink", "purple", "yellow","green"),
border = "black") 
Observam que efectivament, així com es comenta en el treball, el període de confinament va ser de principi de març, on comença a decréixer el valor dels preus, fins a finals de juny, on pareix que s’han tornat a recuperar els valors anteriors a la COVID.
Dividirem la sèrie en tres trams:
Des de 10-2015 fins 12-2018 (serie1_ef), amb el que predirem un període. És a dir predirem de 1-2019 a 12-2019. (amb aquest comprovam que els mètodes funcionen)
Des de 10-2015 fins 12-2019 (serie2_ef), amb la que predirem el període de la COVID (és a dir, l’any 2020)
Des de 10-2015 fins 12-2020 (serie3_ef), amb la que predirem el períodes després de la COVID (és a dir, l’any 2021)
Visualitzem-los:
serie1_ef <- ts(ef$y[1:39],start=c(2015,10),frequency = 12)
serie2_ef <- ts(ef$y[1:51],start=c(2015,10),frequency = 12)
serie3_ef <- ts(ef$y[1:63],start=c(2015,10),frequency = 12)
par(mfrow=c(2,2), mar=c(4,4,4,1)+.1)
plot.ts(serie_ef, main="Serie completa", xlab="Any", ylab="Euros")
plot.ts(serie1_ef, main="Serie abans de la COVID menys un cicle", xlab="Any", ylab="Euros")
plot.ts(serie2_ef, main="Serie abans de la COVID", xlab="Any", ylab="Euros")
plot.ts(serie3_ef, main="Serie abans i durant de la COVID", xlab="Any", ylab="Euros")
2. AJUST DELS MODELS
Abans d’aplicar algun model, estudiem un poc la sèrie:
serie_ef.lm <- lm(y~x, data=ef)
plot.ts(ef$y, main = "Eivissa i Formentera", ylab = "Despeses mensuals en €", xlab="Índex de cada mes")
#dibuixam la recta de regressió sobre els punts
abline(serie_ef.lm, col='red')
Si dibuixam la recta de regressió sobre les nostres dades, tot i que aquestes estan molt disperses i no s’ajusten bé a la recta, podem observar que la tendència és una mica creixent, encara que és manté més o menys constant.
boxplot(serie_ef~cycle(serie_ef), xlab = "mesos", ylab = "Despeses en €", main="Boxplot d'Eivissa i Formentera")
Podem observar també la presència d’estacionalitat, que prenen els seus valors màxims a la temporada d’estiu i els seus mínims en l’hivern, fet que corresponen amb les dades turístiques a les illes.
Aplicarem diversos models per ajustar la nostra sèrie i fer-ne una predicció per llavors comparar quin és el millor.
Vegem com actua cada model:
2.1. ABANSE DE LA COVID
- REGRESSIÓ SEGMENTADA
Com hem comentat anteriorment la recta de regressió no s’ajusta bé a les dades, de fet el valor del \(R^2\) és molt baix:
serie1_ef.lm <- lm(y~x, data=ef[1:39,])
summary(serie1_ef.lm)$adj.r.squared## [1] 0.02707681Per això utilitzam el paquet segmented per ajustar a una
recta de regressió a trossos.
Anem a utilitzar la funció selgmented() per veure quants
de punts de canvi detecta:
punts_canvi_serie1_ef <-selgmented(serie1_ef.lm, Kmax = 10, type="bic")## No. of breakpoints: 2 .. 3 .. 4 .. 5 .. 6 .. 7 .. 8 .. 9 .. 10 ..
##
## BIC to detect no. of breakpoints:
## 0 1 2 3 4 5 6 7
## 419.0407 424.5861 429.3172 432.0721 430.7005 425.9556 386.0211 374.3229
## 8 <NA> <NA>
## 373.3312 374.3312 375.3312
##
## No. of selected breakpoints: 6 (2 breakpoint(s) removed due to small slope diff)Aplicam la funció segmented() que ens calcula la
regressió segmentada:
serie1_ef.seg <- segmented(serie1_ef.lm, seg.Z=~x, psi = c(5,8,16,21,28,34))Aquesta funció ens demana que introduïm els valors on es troben els punts de canvi, i aquesta ens dona el valor estimat. Aquests punts de canvi son:
summary(serie1_ef.seg)$psi## Initial Est. St.Err
## psi1.x 5 5.499227 0.4389605
## psi2.x 8 10.346005 0.4272700
## psi3.x 16 16.491929 0.4791051
## psi4.x 21 23.272538 0.4139897
## psi5.x 28 28.504180 0.3864705
## psi6.x 34 35.159760 0.3795130Que corresponen, aproximadament, a: 2-2016, 7-2016, 1-2017, 8-2017, 2-2018 i 8-2018.
Obtenim que el valor de \(R^2\) per a la regressió segmentada és bastant alt
summary(serie1_ef.seg)$adj.r.squared## [1] 0.8569661Anem a visualitzar la regressió segmentada sobre les nostres dades
#per graficar-ho
p_serie1_ef <- ggplot(ef[1:39,], aes(x=x, y=y)) + geom_line()+
ylim(c(0,1300)) +
ggtitle("Regressió lineal i segmentada sobre les dades \nd'Eivissa i Formentera abans de la COVID") +
xlab('Índex del mes') +
ylab('Despeses mensuals en €')
my.coef1_ef <- coef(serie1_ef.lm) #coeficients de la recta de regressió lineal
p_serie1_ef <- p_serie1_ef + geom_abline(intercept=my.coef1_ef[1], slope=my.coef1_ef[2], color="green") #afegim la recta de regressió lineal
my.model1_ef <- data.frame(x=1:39, y=fitted(serie1_ef.seg)) #model nou amb els valors ajustats de la regressió segmentada
p_serie1_ef <- p_serie1_ef + geom_line(data=my.model1_ef, aes(x=x,y=y), colour="red") #afegim la recta de regressió segmentada
my.lines1_ef <- serie1_ef.seg$psi[,2]# guardam on estan els punts de canvi estimats
p_serie1_ef <- p_serie1_ef + geom_vline(xintercept=my.lines1_ef, linetype="dashed")
p_serie1_ef <- p_serie1_ef + theme(panel.background = element_blank())+ #eliminam el fons i la quadrícula del gràfic
geom_vline(xintercept=0) + geom_hline(yintercept=0)
ggplotly(p_serie1_ef)Visualment també es pot observar que la recta de regressió segmentada s’ajusta millor a les nostres dades.
Anem a calcular les equacions d’aquestes rectes, sabem que les rectes tenen la forma \(y=mx+n\), on \(m\) és la pendent i \(n\) el valor de tall de l’eix y.
Hi ha una funció del paquet segmented que ens dona
aquesta valors de \(n\):
intercept(serie1_ef.seg)## $x
## Est.
## intercept1 933.81
## intercept2 -308.13
## intercept3 1969.50
## intercept4 -1013.20
## intercept5 4245.60
## intercept6 -2504.80
## intercept7 7206.40També tenim una altra funció que ens calcula les pendents:
slope(serie1_ef.seg)## $x
## Est. St.Err. t value CI(95%).l CI(95%).u
## slope1 -92.588 24.405 -3.7938 -142.850 -42.325
## slope2 133.250 24.405 5.4599 82.988 183.510
## slope3 -86.896 18.449 -4.7102 -124.890 -48.901
## slope4 93.960 14.585 6.4423 63.922 124.000
## slope5 -132.000 24.405 -5.4088 -182.270 -81.741
## slope6 104.820 14.585 7.1867 74.779 134.860
## slope7 -171.390 34.514 -4.9656 -242.470 -100.300Aleshores, la nostra recta segmentada és:
\(\hat{y}= \left\{ \begin{array}{lcc} -92.588x + 933.81 & si & x \leq \psi_1 \\ \\ 133.250x - 308.13 & si & \psi_1 < x \leq \psi_2 \\ \\ -86.896x + 1969.50 & si & \psi_2 < x \leq \psi_3 \\ \\ 93.960x - 1013.20 & si & \psi_3 < x \leq \psi_4 \\ \\ -132.000x + 4245.60 & si & \psi_4 < x \leq \psi_5 \\ \\ 104.820x - 2504.80 & si & \psi_5 < x \leq \psi_6 \\ \\ -171.390x + 7206.40 & si & \psi_6 < x \\ \end{array} \right.\)
on \(\psi_1= 5.5\), \(\psi_2= 10.35\), \(\psi_3= 16.49\), \(\psi_4=23.27\), \(\psi_5= 28.5\) i \(\psi_6= 35.16\).
Vegem els errors del model (ens interessa el RMSE):
accuracy(serie1_ef.seg)## ME RMSE MAE MPE MAPE MASE
## Training set 4.372882e-15 61.79027 46.87676 -0.612519 6.148984 0.2692355Anem a fer la predicció de 1-2019 a 12-2019.
El darrer punt de canvi que tenim és en agost de 2018 i, seguint el mateix criteri que amb les dades de Mallorca, podem considerar que els següents es donaran en 2-2019, 8-2019 i 2-2020. Necessitam calcular les pendents de les rectes d’entre agost de 2018 i febrer de 2020, per poder fer la predicció i calcular l’error.
Per calcular les pendents de les noves rectes farem la mitjana de les pendents anteriors. De les pendents ja calculades en el model obviarem la primera i la darrera, ja que no són vàlides. Per tant,
La pendent de 8-2018 a 2-2019 i de 8-2019 a 2-2020 serà : -109.448
La pendent de 2-2019 a 8-2019 serà : 110.677
Ara, seguint el mateix procediment que en el cas de Mallorca, els nous punts d’intersecció són: 5028.86, -3996.26 i 6349.61.
Llavors l’equació per a la predicció és:
\(\hat{y}= \left\{ \begin{array}{lcc} -109.448x + 5028.86 & si & x \leq \psi_7 \\ \\ 110.677x - 3996.26 & si & \psi_7 < x \leq \psi_8 \\ \\ -109.448x + 6349.61 & si & \psi_8 < x \\ \end{array} \right.\)
on \(\psi_7 = 41\) i \(\psi_8 = 47\).
# per graficar-la
p1_serie_ef <- ggplot(ef, aes(x=x, y=y)) + geom_line()+
ylim(c(0,1300)) +
ggtitle("Predicció d'e les Illes Balears'Eivissa i Formentera amb el \nmodel de regressió segmentada abans de la COVID") +
xlab('Índex del mes') +
ylab('Despeses mensuals en €')
my.model1_ef <- data.frame(x=1:39, y=fitted(serie1_ef.seg)) #model nou amb els valors ajustats de la regressió segmentada
p1_serie_ef <- p1_serie_ef + geom_line(data=my.model1_ef, aes(x=x,y=y), colour="red") #afegim la recta de regressió segmentada
my.lines1_ef <- serie1_ef.seg$psi[,2]# guardam on estan els punts de canvi estimats
p1_serie_ef <- p1_serie_ef + geom_vline(xintercept=my.lines1_ef, linetype="dashed")
p1_serie_ef <- p1_serie_ef + geom_vline(xintercept=0) + geom_hline(yintercept=0) + theme(panel.background = element_blank()) #eliminam el fons i la quadrícula del gràfic
p1_serie_ef <- p1_serie_ef + geom_abline(intercept = 5028.86, slope=-109.448, colour="green") +
geom_abline(intercept = -3996.26, slope=110.677, colour="blue") +
geom_abline(intercept = 6349.61, slope=-109.448, colour="orange")
ggplotly(p1_serie_ef)Calculem l’error de la predicció:
o1_ef<-c(serie_ef[40:51]) #observacions reals de la sèrie de gener de 2019 a desembre de 2019
f0_ef <- -109.448 * 40 + 5028.68 #predicció de gener 2019
v1_ef=c(41:47)
f1_ef <- sapply(v1_ef, function(x) 110.677*x-3996.26) #predicció de febrer 2019 a agost 2019
v2_ef=c(48:51)
f2_ef <- sapply(v2_ef, function(x) -109.448*x + 6349.61 )#predicció de agost 2019 a desembre 2019
p1_ef<- c(f0_ef,f1_ef,f2_ef) #predicció de 2019
sqrt(sum((p1_ef-o1_ef)^2)/12) #error de la predicció## [1] 109.3058- DECOMPOSE
Recordem la sèrie
plot.ts(serie1_ef, main= "Eivissa i Formentera abans de la COVID", xlab="Any", ylab="Despeses mensuals en €")
Ja hem dit anteriorment que té una tendència mínima i podem observar també que no hi ha variabilitat.
El mètode clàssic de descomposició separa la sèrie en subseries corresponents a la tendència, la estacionalitat i el renou.
Aquestes components es poden combinar de manera additiva o multiplicativa. En el nostre cas utilitzam el model additiu: \(y_t = \mu_t + S_t + a_t\)
decompose_s1_ef<-decompose(serie1_ef)
plot(decompose_s1_ef, xlab="Any")
El decompose, per calcular aquestes noves tendències, el
que fa és agafar les sis tendències anteriors i les sis posteriors de la
sèrie original i en fa la mitjana. És per això que la primera que
obtenim és en abril de 2016 i la darrera en juny de 2018. Notem que per
a la predicció ens quedarem amb el valor de la darrera tendència del
decompose.
t1_ef <- decompose_s1_ef$trend #tendències de la sèrie sense el tros a predir abans de la COVID
t1_ef## Jan Feb Mar Apr May Jun Jul Aug
## 2015
## 2016 NA NA NA 743.1446 747.8983 758.9400 759.6646 763.9542
## 2017 806.7158 817.0208 824.6025 837.5362 845.8621 843.2208 842.9546 837.8417
## 2018 831.0383 830.7542 836.4517 830.5658 816.6008 804.8000 NA NA
## Sep Oct Nov Dec
## 2015 NA NA NA
## 2016 773.8896 782.0333 791.3829 797.4979
## 2017 833.9504 838.2292 840.2021 837.6517
## 2018 NA NA NA NAEls valors de les components estacionals els calcula fent la mitjana per mesos, és a dir, per calcular la componen de gener, agafa els valors de tots els geners que tenim a la sèrie original i en fa la mitjana. Per tant, només tenim 12 valors, un per a cada mes.
s1_ef <- decompose_s1_ef$seasonal
s1_ef <- s1_ef[4:15] #estacionalitat de gener a desembre
s1_ef## [1] -306.83481 -201.66023 -246.34481 -95.47495 19.08352 117.80366
## [7] 299.75769 236.97936 121.06727 197.11102 -22.20023 -119.28752Anem a fer la predicció d’aquest model:
a1_ef <- c(s1_ef[7:12],s1_ef) #estacionalitat de juliol 2018 a desembre 2019 (es per poder fer la predicció)
pred1_decompose_ef <- sapply(a1_ef, function(x) 804.8000 + x) #predicció de juliol de 2018 a desembre 2019p1_dec_ef<-c(serie1_ef[1:33], pred1_decompose_ef)
A1_ef<- data.frame("x" = ef[1:51,]$x, "y" = ef[1:51,]$y, "p"= p1_dec_ef)
grafica1_ef_dec <- ggplot(A1_ef)+
geom_line(aes(x,p), color="red")+
geom_line(aes(x,y))+
geom_vline(xintercept=0) + geom_hline(yintercept=0)+
labs(title="Predicció d'Eivissa i Formntera abans de la COVID amb el model \nde descomposició", x="Índex del mes", y="Despeses mensuals en €")+ theme(panel.background = element_blank())
grafica1_ef_dec
Podem observar que aquest model fa una predicció bastant bona. Calculem l’error de la predicció:
(Notem que la predicció és des de juliol de 2018 a desembre de 2019 i per calcular l’error només volem de gener de 2019 a desembre de 2019.)
sqrt(sum((c(serie_ef[40:51]- pred1_decompose_ef[7:18]))^2)/12) #error predicció de gener a desembre de 2019## [1] 84.1187Notem que també tenim una altra instrucció a R per fer prediccions
d’una sèrie, la funció predict(). Aquesta és basa en un
model ETS, anem a veure la predicció:
prediccio1_ef <- predict(serie1_ef,h=12)
plot(prediccio1_ef, xlab="Any", ylab="Despeses mensuals en €")
Pareix que la predicció és bastant bona, ja que el cicle predit segueix un mateix patró que els anteriors. Anem a veure la predicció juntament amb la sèrie original:
df_prediccio1_ef <- data.frame(prediccio1_ef)
p1_ets_ef <- data.frame("x"= 1:51, "PointForecast"=c(serie1_ef,df_prediccio1_ef$Point.Forecast), "Lo80" = c(rep(NA,39),df_prediccio1_ef$Lo.80), "Hi80" = c(rep(NA,39),df_prediccio1_ef$Hi.80), "Lo95" = c(rep(NA,39),df_prediccio1_ef$Lo.95),"Hi95" = c(rep(NA,39),df_prediccio1_ef$Hi.95))
grafica_pred1_ets <- ggplot((ef[1:51,]))+
geom_ribbon(data = p1_ets_ef, aes(x, ymin = Lo95, ymax = Hi95), fill = "blue", alpha = 0.25) +
geom_ribbon(data =p1_ets_ef, aes(x, ymin = Lo80, ymax = Hi80), fill = "blue", alpha = 0.25) +
geom_line(data = p1_ets_ef, aes(x, PointForecast), colour = "blue") +
geom_line(aes(x,y), color="red")+
labs(title="Predicció del model ETS(M,N,M) a Eivissa i Formentera \nabans de la COVID", x="Índex del mes", y="Despeses mensuals en €") +
geom_vline(xintercept=0) + geom_hline(yintercept=0)+ theme(panel.background = element_blank())
grafica_pred1_ets
Podem observar que la predicció és bastant bona, ja que continua seguint un mateix patró. I l’error comés és d’uns 76 euros.
accuracy(prediccio1_ef,serie2_ef)## ME RMSE MAE MPE MAPE MASE
## Training set -0.9578372 60.39005 47.93061 -0.9764864 6.671517 0.4957655
## Test set 30.7451149 76.44673 66.73647 1.9682600 8.708841 0.6902821
## ACF1 Theil's U
## Training set 0.1630727 NA
## Test set 0.2101528 0.4449581- SARIMA
Anem a veure quin model obtenim considerant un model estacional
Pel que hem vist anteriorment, podem considerar que no hi ha tendència, llavors no fa falta fer cap diferència a la part regular, no obstant, sí que cal fer una diferenciació d’orde 12. Vegem l’ACF i el PACF de la sèrie a predir:
par(mfrow=c(1,2))
acf(serie1_ef)
pacf(serie1_ef)
Per a la part regular obtenim un ARIMA(1,0,2) Feim una diferenciació a la part estacional, és a dir, d’ordre 12
serie1_ef_dif <- diff(serie1_ef,12)
plot(serie1_ef_dif, main="Sèrie sense estacionalitat", xlab="Any", ylab="Sèrie diferenciada")
Aquesta és la sèrie sense estacionalitat ni tendència, vegem com es modifiquen l’ACF i el PACF.
par(mfrow=c(2,1))
acf(serie1_ef_dif, lag=36)
pacf(serie1_ef_dif,lag=36)
Per a la part estacional obtenim que P=0, D=1, Q=0
És a dir, obtenim un ARIMA(1,0,2)(0,1,0)[12], vegem quin model els proposa R:
auto.arima(serie1_ef)## Series: serie1_ef
## ARIMA(0,1,1)(1,1,0)[12]
##
## Coefficients:
## ma1 sar1
## -0.6686 -0.4226
## s.e. 0.1826 0.2143
##
## sigma^2 = 9792: log likelihood = -156.79
## AIC=319.59 AICc=320.68 BIC=323.36R ens proposa un ARIMA(0,1,1)(1,1,0)[12]. Per tant les propostes de model ARIMA són:
model1_ef<-arima(serie1_ef, order=c(1,0,2), seasonal = list(order=c(0,1,0), period=12))
model2_ef <- arima(serie1_ef, order=c(0,1,1), seasonal = list( order=c(1,1,0), period=12))Amb uns errors d’uns 85 i 78 euros cada un.
accuracy(model1_ef)## ME RMSE MAE MPE MAPE MASE
## Training set -1.700455 84.64409 59.40385 -1.134031 7.905317 0.4795375
## ACF1
## Training set 0.04145267accuracy(model2_ef)## ME RMSE MAE MPE MAPE MASE ACF1
## Training set -9.59489 77.62801 54.45449 -2.184592 7.553735 0.4395838 0.08070217Anem a fer la predicció de 2019 d’aquests models:
fc_m1_ef <-forecast(model1_ef, h=12)
fc_m2_ef <- forecast(model2_ef,h=12)Visualitzem-les
#Model 1
pre1_ef <- data.frame("Point Forecast" = serie1_ef, "Lo 80" = rep(NA,39), "Hi 80"= rep(NA,39), "Lo 95" = rep(NA,39), "Hi 95" = rep(NA,39)) #dades abans de la predicció amb NA als intervals ja que només els volem per la predicció
pred_m1_ef <-data.frame(fc_m1_ef)
grafica_m1_ef <- data.frame("x" = 1:51, "PointForecast" = c(pre1_ef$Point.Forecast,pred_m1_ef$Point.Forecast), "Lo80" = c(pre1_ef$Lo.80, pred_m1_ef$Lo.80), "Hi80"= c(pre1_ef$Hi.80, pred_m1_ef$Hi.80), "Lo95" = c(pre1_ef$Lo.95, pred_m1_ef$Lo.95), "Hi95" = c(pre1_ef$Hi.95, pred_m1_ef$Hi.95))
grafica1_ef <- ggplot(ef[1:51,])+
geom_ribbon(data = grafica_m1_ef, aes(x, ymin = Lo95, ymax = Hi95), fill = "blue", alpha = 0.25) +
geom_ribbon(data = grafica_m1_ef, aes(x, ymin = Lo80, ymax = Hi80), fill = "blue", alpha = 0.25) +
geom_line(data = grafica_m1_ef, aes(x, PointForecast), colour = "blue") +
geom_line(aes(x,y), color="red")+
labs(title="Predicció d'Eivissa i Formentera amb el model ARIMA(1,0,2)(0,1,0)[12] \nabans de la COVID", x="Índex del mes", y="Despeses mensuals en €") +
geom_vline(xintercept=0) + geom_hline(yintercept=0)+ theme(panel.background = element_blank())
grafica1_ef
#Model 2
pred_m2_ef <-data.frame(fc_m2_ef)
grafica_m2_ef <- data.frame("x" = 1:51, "PointForecast" = c(pre1_ef$Point.Forecast,pred_m2_ef$Point.Forecast), "Lo80" = c(pre1_ef$Lo.80, pred_m2_ef$Lo.80), "Hi80"= c(pre1_ef$Hi.80, pred_m2_ef$Hi.80), "Lo95" = c(pre1_ef$Lo.95, pred_m2_ef$Lo.95), "Hi95" = c(pre1_ef$Hi.95, pred_m2_ef$Hi.95))
grafica2_ef <- ggplot(ef[1:51,])+
geom_ribbon(data = grafica_m2_ef, aes(x, ymin = Lo95, ymax = Hi95), fill = "blue", alpha = 0.25) +
geom_ribbon(data = grafica_m2_ef, aes(x, ymin = Lo80, ymax = Hi80), fill = "blue", alpha = 0.25) +
geom_line(data = grafica_m2_ef, aes(x, PointForecast), colour = "blue") +
geom_line(aes(x,y), color="red")+
labs(title="Predicció d'Eivissa i Formentera amb el model ARIMA (0,1,1)(1,1,0)[12] \nabans de la COVID", x="Índex del mes", y="Despeses mensuals en €") +
geom_vline(xintercept=0) + geom_hline(yintercept=0)+ theme(panel.background = element_blank())
grafica2_ef
Vegem i estudiem els errors de la predicció:
accuracy(fc_m1_ef, serie_ef[40:51], h=12)## ME RMSE MAE MPE MAPE MASE
## Training set -1.700455 84.64409 59.40385 -1.134031 7.905317 0.4795375
## Test set 52.755681 80.73996 64.25420 6.415168 8.130005 0.5186920
## ACF1
## Training set 0.04145267
## Test set NAaccuracy(fc_m2_ef,serie_ef[40:51], h=12)## ME RMSE MAE MPE MAPE MASE ACF1
## Training set -9.59489 77.62801 54.45449 -2.184592 7.553735 0.4395838 0.08070217
## Test set 68.97582 93.75512 78.51529 8.037159 9.918771 0.6338147 NAcheckresiduals(fc_m1_ef)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(1,0,2)(0,1,0)[12]
## Q* = 5.3679, df = 5, p-value = 0.3727
##
## Model df: 3. Total lags used: 8checkresiduals(fc_m2_ef)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(0,1,1)(1,1,0)[12]
## Q* = 6.655, df = 6, p-value = 0.3539
##
## Model df: 2. Total lags used: 8par(mfrow=c(1,2))
qqPlot(fc_m1_ef$residuals, id=FALSE, mean=mean(fc_m1_ef$residuals), sd=sd(fc_m1_ef$residuals))
qqPlot(fc_m2_ef$residuals, id=FALSE, mean=mean(fc_m2_ef$residuals), sd=sd(fc_m2_ef$residuals))
shapiro.test(fc_m1_ef$residuals)##
## Shapiro-Wilk normality test
##
## data: fc_m1_ef$residuals
## W = 0.93959, p-value = 0.03689shapiro.test(fc_m2_ef$residuals)##
## Shapiro-Wilk normality test
##
## data: fc_m2_ef$residuals
## W = 0.94603, p-value = 0.06045lillie.test(fc_m1_ef$residuals)##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: fc_m1_ef$residuals
## D = 0.2282, p-value = 2.121e-05lillie.test(fc_m2_ef$residuals)##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: fc_m2_ef$residuals
## D = 0.19193, p-value = 0.0009128Resum dels errors que cometen cada un dels models anteriors:
| reg. segmentada | descomposició clàssica | ETS(M,N,M) | ARIMA (1,0,2)(0,1,0)[12] | ARIMA (0,1,1)(1,1,0)[12] | |
|---|---|---|---|---|---|
| error model | 61.79 | 60.39 | 84.64 | 77.62 | |
| error predicció | 109.3058 | 84.1187 | 76.45 | 80.74 | 93.76 |
2.2. DURANT LA COVID
- REGRESSIÓ SEGMENTADA
El valor de \(R^2\) per a la regressió lineal és:
serie2_ef.lm <- lm(y~x, data= ef[1:51,])
summary(serie2_ef.lm)$adj.r.squared## [1] 0.007558416Aleshores utilitzem el paquet segmented per ajustar les
nostres dades a una regressió lineal segmentada i millorar aquest valor.
Vegem els punts de canvi:
punts_canvi_serie2_ef <-selgmented(serie2_ef.lm, Kmax = 10, type="bic")## No. of breakpoints: 2 .. 3 .. 4 .. 5 .. 6 .. 7 .. 8 .. 9 .. 10 ..
##
## BIC to detect no. of breakpoints:
## 0 1 2 3 4 5 6 7
## 551.5657 557.6175 564.7267 572.0761 573.9956 563.8713 564.6504 530.2193
## 8 9 10
## 508.0282 504.6689 508.9345
##
## No. of selected breakpoints: 8 (1 breakpoint(s) removed due to small slope diff)serie2_ef.seg <- segmented(serie2_ef.lm, seg.Z=~x, psi = c(5,10,14,21,29,35,40,47))
summary(serie2_ef.seg)$psi## Initial Est. St.Err
## psi1.x 5 5.499227 0.4621615
## psi2.x 10 10.346005 0.4498531
## psi3.x 14 16.491929 0.5044278
## psi4.x 21 23.272538 0.4358708
## psi5.x 29 28.484784 0.4244957
## psi6.x 35 34.969280 0.4036198
## psi7.x 40 40.649312 0.3482145
## psi8.x 47 46.591245 0.3590303Aquests són en: 3-2016, 7-2016, 1-2017, 8-2017, 1-2018, 8-2018, 2-2019 i 8-2019.
Ara el valor de \(R^2\) de la segmentada és
summary(serie2_ef.seg)$adj.r.squared## [1] 0.852397Anem a visualitzar la regressió segmentada sobre les nostres dades
p_serie2_ef <- ggplot(ef[1:51,], aes(x=x, y=y)) + geom_line()+
ylim(c(0,1300))+
ggtitle("Regressió lineal i segmentada sobre les dades \nd'Eivissa i Formentera durant la COVID") + xlab('índex del mes')+ ylab("Despeses mensuals en €")
my.coef2_ef <- coef(serie2_ef.lm) #coeficients de la recta de regressió lineal
p_serie2_ef <- p_serie2_ef + geom_abline(intercept=my.coef2_ef[1], slope=my.coef2_ef[2], color="green") #afegim la recta de regressió lineal
my.model2_ef <- data.frame(x=1:51, y=fitted(serie2_ef.seg)) #model nou amb els valors ajustats de la regressió segmentada
p_serie2_ef <- p_serie2_ef + geom_line(data=my.model2_ef, aes(x=x,y=y), colour="red") #afegim la recta de regressió segmentada
my.lines2_ef <- serie2_ef.seg$psi[,2]# guardam on estan els punts de canvi estimats
p_serie2_ef <- p_serie2_ef + geom_vline(xintercept=my.lines2_ef, linetype="dashed")
p_serie2_ef <- p_serie2_ef + geom_vline(xintercept=0) + geom_hline(yintercept=0) + theme(panel.background = element_blank()) #eliminam el fons i la quadrícula del gràfic
ggplotly(p_serie2_ef)Per poder escriure la funció necessitam les pendents i els punts d’intersecció amb l’eix OY, que ens ho donen les següents funcions:
#punts d'intersecció
intercept(serie2_ef.seg)## $x
## Est.
## intercept1 933.81
## intercept2 -308.13
## intercept3 1969.50
## intercept4 -1013.20
## intercept5 4245.60
## intercept6 -2440.20
## intercept7 5850.90
## intercept8 -5229.00
## intercept9 7654.80#pendents
slope(serie2_ef.seg)## $x
## Est. St.Err. t value CI(95%).l CI(95%).u
## slope1 -92.588 25.695 -3.6033 -144.870 -40.311
## slope2 133.250 25.695 5.1858 80.974 185.530
## slope3 -86.896 19.424 -4.4737 -126.410 -47.378
## slope4 93.960 15.356 6.1189 62.719 125.200
## slope5 -132.000 25.695 -5.1373 -184.280 -79.727
## slope6 102.710 19.424 5.2880 63.194 142.230
## slope7 -134.390 19.424 -6.9187 -173.900 -94.868
## slope8 138.190 19.424 7.1144 98.670 177.710
## slope9 -138.340 25.695 -5.3839 -190.620 -86.063L’error del model de regressió segmentada és (ens interessa el RMSE):
accuracy(serie2_ef.seg)## ME RMSE MAE MPE MAPE MASE
## Training set -3.342914e-15 65.36159 51.33711 -0.6744592 6.975136 0.2810544Anem a fer la predicció, recordem que els punts de canvi són en 3-2016, 7-2016, 1-2017, 8-2017, 1-2018, 8-2018, 2-2019 i 8-2019. Aleshores, igual que abans els següents punts de canvi es donen en gener, agost i gener, per tant:
- La pendent de 8-2019 a 1-2020 i de 8-2020 a 1-2021 serà: -117.762
- La pendent de 1-2020 a 8-2020 serà : 117.0275
Seguin el mateix procediment que amb les dades de Mallorca, els nous punts d’intersecció són 6695.8, -5513.25 i 8339.33.
És a dir, la funció per a la predicció de 2020 és:
\(\hat{y}= \left\{ \begin{array}{lcc} -117.762x + 6695.8 & si & x \leq \psi_9 \\ \\ 117.0275x - 5513.25 & si & \psi_9 < x \leq \psi_{10} \\ \\ -117.762x + 8339.33& si & \psi_{10} < x \\ \end{array} \right.\)
on \(\psi_7 = 52\) i \(\psi_8 = 59\).
#ho graficam
p2_serie_ef <- ggplot(ef, aes(x=x, y=y)) + geom_line()+
ylim(c(0,1400))+
ggtitle("Predicció d'Eivissa i Formentera amb el model \nde regressió segmentada durant la COVID") +
xlab('Índex del mes') +
ylab('Despeses mensuals en €')
my.model2_ef <- data.frame(x=1:51, y=fitted(serie2_ef.seg)) #model nou amb els valors ajustats de la regressiño segmentada
p2_serie_ef <- p2_serie_ef + geom_line(data=my.model2_ef, aes(x=x,y=y), colour="red") #afegim la recta de regressió segmentada
my.lines2_ef <- serie2_ef.seg$psi[,2]# guardam on estan els punts de canvi estimats
p2_serie_ef <- p2_serie_ef + geom_vline(xintercept=my.lines2_ef, linetype="dashed")
p2_serie_ef <- p2_serie_ef + theme(panel.background = element_blank()) + #eliminam el fons i la quadrícula del gràfic
geom_vline(xintercept=0) + geom_hline(yintercept=0)
p2_serie_ef <- p2_serie_ef + geom_abline(intercept = 6695.8, slope=-117.762, colour="green") +
geom_abline(intercept = -5513.25 , slope=117.0275, colour="blue") +
geom_abline(intercept =8339.33 , slope=-117.762, colour="orange")
ggplotly(p2_serie_ef)Calculem l’error de la predicció:
o2_ef<-c(serie_ef[52:63]) # dades reals per fer predicció del 2020
v3_ef=c(52:59)
f3_ef <- sapply(v3_ef, function(x) 117.0275*x-5513.25) #predicció de gener 2020 a agost 2020
v4_ef=c(60:63)
f4_ef <- sapply(v4_ef, function(x) -117.762*x + 8339.33) #predicció de setembre 2020 a desembre 2020
p2_ef <- c(f3_ef,f4_ef) #predicció de 2020
sqrt(sum((p2_ef-o2_ef)^2)/12)## [1] 530.5246- DECOMPOSE
La descomposició de la sèrie en aquest cas és la següent
decompose_s2_ef <- decompose(serie2_ef)
plot(decompose_s2_ef, xlab="Any")
On les components de tendència són:
t2_ef <- decompose_s2_ef$trend #tendències de la sèrie sense el tros a predir abans de la COVID
t2_ef## Jan Feb Mar Apr May Jun Jul Aug
## 2015
## 2016 NA NA NA 743.1446 747.8983 758.9400 759.6646 763.9542
## 2017 806.7158 817.0208 824.6025 837.5362 845.8621 843.2208 842.9546 837.8417
## 2018 831.0383 830.7542 836.4517 830.5658 816.6008 804.8000 801.1833 800.7592
## 2019 793.0312 797.2383 795.3000 794.5121 794.6762 798.6983 NA NA
## Sep Oct Nov Dec
## 2015 NA NA NA
## 2016 773.8896 782.0333 791.3829 797.4979
## 2017 833.9504 838.2292 840.2021 837.6517
## 2018 794.4083 787.4583 785.2383 787.7225
## 2019 NA NA NA NAI les estacionals:
s2_ef <- decompose_s2_ef$seasonal #estacionalitat
s2_ef <- s2_ef[4:15] #estacionalitat de gener a desembreFem la predicció:
a2_ef <- c(s2_ef[7:12],s2_ef) #estacionalitat de juliol 2019 a desembre 2029 (es per poder fer la predicció)
pred2_decompose_ef <- sapply(a2_ef, function(x) 798.6983 + x) #predicció de juliol de 2019 a desembre 2020
p2_dec_ef<-c(serie2_ef[1:45], pred2_decompose_ef)
A2_ef<- data.frame("x" = ef[1:63,]$x, "y" = ef[1:63,]$y, "p"= p2_dec_ef)
grafica2_ef_dec <- ggplot(A2_ef)+
geom_line(aes(x,p), color="red")+
geom_line(aes(x,y))+
labs(title="Predicció d'Eivissa i Formentera durant la COVID amb el model \nde descomposició", x="Índex del mes", y="Despesese mensuals en €")+
geom_vline(xintercept = 0)+ geom_hline(yintercept = 0)+ theme(panel.background = element_blank())
grafica2_ef_dec
L’error que es comet és:
sqrt(sum((c(serie_ef[52:63]- pred2_decompose_ef[7:18]))^2)/12)## [1] 348.4049Vegem, igual que abans, la predicció amb la funció
predict():
prediccio2_ef <- predict(serie2_ef,h=12)
plot(prediccio2_ef, xlab="Any", ylab ="Despeses menusals en €")
Vegem com queda la predicció sobre la sèrie original:
df_prediccio2_ef <- data.frame(prediccio2_ef)
p2_ets_ef <- data.frame("x"= 1:63, "PointForecast"=c(serie2_ef,df_prediccio2_ef$Point.Forecast), "Lo80" = c(rep(NA,51),df_prediccio2_ef$Lo.80), "Hi80" = c(rep(NA,51),df_prediccio2_ef$Hi.80), "Lo95" = c(rep(NA,51),df_prediccio2_ef$Lo.95),"Hi95" = c(rep(NA,51),df_prediccio2_ef$Hi.95))
grafica_pred2_ets <- ggplot((ef[1:63,]))+
geom_ribbon(data = p2_ets_ef, aes(x, ymin = Lo95, ymax = Hi95), fill = "blue", alpha = 0.25) +
geom_ribbon(data =p2_ets_ef, aes(x, ymin = Lo80, ymax = Hi80), fill = "blue", alpha = 0.25) +
geom_line(data = p2_ets_ef, aes(x, PointForecast), colour = "blue") +
geom_line(aes(x,y), color="red")+
labs(title="Predicció del model ETS(A,N,A) a Eivissa i Formentera durant la COVID", y="Despeses mensuals en €", x="Índex del mes") +
geom_vline(xintercept = 0) + geom_hline(yintercept = 0)+ theme(panel.background = element_blank())
grafica_pred2_ets
Si calculam l’error comés, aquest és d’uns 350 euros.
accuracy(prediccio2_ef, serie3_ef)## ME RMSE MAE MPE MAPE MASE
## Training set 1.253437 63.58661 49.70861 -0.6009612 6.749908 0.5855261
## Test set -230.523743 350.31320 232.13630 -Inf Inf 2.7343725
## ACF1 Theil's U
## Training set 0.1148920 NA
## Test set 0.4007914 0- SARIMA
Quin model proposam nosaltres? Vegem l’ACF i el PACF:
par(mfrow=c(1,2))
acf(serie2_ef)
pacf(serie2_ef)
Per a la part regular obtenim que p=1, q=2 i d=0.
Observam que hi ha estacionalitat, llavors feim una diferenciació d’ordre 12.
serie2_ef_diff <- diff(serie2_ef,12)
plot(serie2_ef_diff, main="Sèrie sense estacionalitat", xlab="Any", ylab="Sèrie diferenciada")
Ara ja no s’observa l’estacionalitat, llavors hem de fer una diferenciació D=1. Vegem quins són els nous ACF i PACF.
par(mfrow=c(1,2))
acf(serie2_ef_diff,lag=36)
pacf(serie2_ef_diff,lag=36)
Obtenim que P=0 i Q=0. Per tant, el model que nosaltres proposam es un ARIMA(1,0,2)(0,1,0)[12]
R proposa el següent model:
auto.arima(serie2_ef)## Series: serie2_ef
## ARIMA(0,1,1)(0,1,1)[12]
##
## Coefficients:
## ma1 sma1
## -0.7560 -0.5981
## s.e. 0.1256 0.3089
##
## sigma^2 = 7322: log likelihood = -224.98
## AIC=455.95 AICc=456.66 BIC=460.87Els dos models que tenim són
model3_ef <- arima(serie2_ef, order=c(1,0,2), seasonal = list(order=c(0,1,0), period=12))
model4_ef <- arima(serie2_ef, order=c(0,1,1), seasonal = list( order=c(0,1,1), period=12))I els seus errors
accuracy(model3_ef)## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 6.692082 82.17198 60.48591 0.0511561 8.013875 0.4862884 0.02451003accuracy(model4_ef)## ME RMSE MAE MPE MAPE MASE
## Training set -5.597112 71.89198 53.10005 -1.708372 7.395149 0.4269083
## ACF1
## Training set 0.08132473Les prediccions són:
fc_m3_ef <- forecast(model3_ef, h=12)
fc_m4_ef <- forecast(model4_ef,h=12)#primer model
pre2_ef <- data.frame("Point Forecast" = serie2_ef, "Lo 80" = rep(NA,51), "Hi 80"= rep(NA,51), "Lo 95" = rep(NA,51), "Hi 95" = rep(NA,51))
pred_m3_ef <-data.frame(fc_m3_ef)
grafica_m3_ef <- data.frame("x" = 1:63, "PointForecast" = c(pre2_ef$Point.Forecast,pred_m3_ef$Point.Forecast), "Lo80" = c(pre2_ef$Lo.80, pred_m3_ef$Lo.80), "Hi80"= c(pre2_ef$Hi.80, pred_m3_ef$Hi.80), "Lo95" = c(pre2_ef$Lo.95, pred_m3_ef$Lo.95), "Hi95" = c(pre2_ef$Hi.95, pred_m3_ef$Hi.95))
grafica3_ef <- ggplot(ef[1:63,])+
geom_ribbon(data = grafica_m3_ef, aes(x, ymin = Lo95, ymax = Hi95), fill = "blue", alpha = 0.25) +
geom_ribbon(data = grafica_m3_ef, aes(x, ymin = Lo80, ymax = Hi80), fill = "blue", alpha = 0.25) +
geom_line(data = grafica_m3_ef, aes(x, PointForecast), colour = "blue") +
geom_line(aes(x,y), color="red")+
labs(title="Prediccó d'Eivissa i Formentera amb el model ARIMA (1,0,2)(0,1,0)[12] \ndurant la COVID", x="Índex del mes", y="Despeses mensuals en €")+
geom_vline(xintercept=0) + geom_hline(yintercept=0) + theme(panel.background = element_blank())
grafica3_ef
#segon model
pred_m4_ef <-data.frame(fc_m4_ef)
grafica_m4_ef <- data.frame("x" = 1:63, "PointForecast" = c(pre2_ef$Point.Forecast,pred_m4_ef$Point.Forecast), "Lo80" = c(pre2_ef$Lo.80, pred_m4_ef$Lo.80), "Hi80"= c(pre2_ef$Hi.80, pred_m4_ef$Hi.80), "Lo95" = c(pre2_ef$Lo.95, pred_m4_ef$Lo.95), "Hi95" = c(pre2_ef$Hi.95, pred_m4_ef$Hi.95))
grafica4_ef <- ggplot(ef[1:63,])+
geom_ribbon(data = grafica_m4_ef, aes(x, ymin = Lo95, ymax = Hi95), fill = "blue", alpha = 0.25) +
geom_ribbon(data = grafica_m4_ef, aes(x, ymin = Lo80, ymax = Hi80), fill = "blue", alpha = 0.25) +
geom_line(data = grafica_m4_ef, aes(x, PointForecast), colour = "blue") +
geom_line(aes(x,y), color="red")+
labs(title="Predicció d'Eivissa i Formentera amb el model ARIMA (0,1,1)(0,1,1)[12] \ndurant la COVID", x="Índex del mes", y="Despeses mensuals en €")+
geom_vline(xintercept=0) + geom_hline(yintercept=0) + theme(panel.background = element_blank())
grafica4_ef
Vegem i estudiem els seus errors:
accuracy(fc_m3_ef,serie_ef[52:63], h=12)## ME RMSE MAE MPE MAPE MASE
## Training set 6.692082 82.17198 60.48591 0.0511561 8.013875 0.4862884
## Test set -248.723052 374.77285 259.75035 -Inf Inf 2.0883141
## ACF1
## Training set 0.02451003
## Test set NAaccuracy(fc_m4_ef,serie_ef[52:63], h=12)## ME RMSE MAE MPE MAPE MASE
## Training set -5.597112 71.89198 53.10005 -1.708372 7.395149 0.4269083
## Test set -241.928797 360.06127 241.92880 -Inf Inf 1.9450342
## ACF1
## Training set 0.08132473
## Test set NAcheckresiduals(fc_m3_ef)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(1,0,2)(0,1,0)[12]
## Q* = 12.22, df = 7, p-value = 0.09355
##
## Model df: 3. Total lags used: 10checkresiduals(fc_m4_ef)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(0,1,1)(0,1,1)[12]
## Q* = 7.1127, df = 8, p-value = 0.5245
##
## Model df: 2. Total lags used: 10par(mfrow=c(1,2))
qqPlot(fc_m3_ef$residuals, id=FALSE, mean=mean(fc_m3_ef$residuals), sd=sd(fc_m3_ef$residuals))
qqPlot(fc_m4_ef$residuals, id=FALSE, mean=mean(fc_m4_ef$residuals), sd=sd(fc_m4_ef$residuals))
shapiro.test(fc_m3_ef$residuals)##
## Shapiro-Wilk normality test
##
## data: fc_m3_ef$residuals
## W = 0.96971, p-value = 0.2151shapiro.test(fc_m4_ef$residuals)##
## Shapiro-Wilk normality test
##
## data: fc_m4_ef$residuals
## W = 0.95558, p-value = 0.0541lillie.test(fc_m3_ef$residuals)##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: fc_m3_ef$residuals
## D = 0.15645, p-value = 0.003194lillie.test(fc_m4_ef$residuals)##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: fc_m4_ef$residuals
## D = 0.1427, p-value = 0.01114Resum dels errors que comet cada model:
| reg. segmentada | descomposició clàssica | ETS(A,N,A) | ARIMA (1,0,2)(0,1,0)[12] | ARIMA (0,1,1)(0,1,1)[12] | |
|---|---|---|---|---|---|
| error model | 65.36 | 63.59 | 82.17 | 71.84 | |
| error predicció | 530.5246 | 348.405 | 350.31 | 374.77 | 360.06 |
2.3. DESPRÉS DE LA COVID
- REGRESSIÓ SEGMENTADA
El valor de \(R^2\) per a la regressió lineal és molt baix
serie3_ef.lm <- lm(y~x, data= ef[1:63,])
summary(serie3_ef.lm)$adj.r.squared## [1] 0.00449053Anem a fer ús del paquet segmented.
Vegem els punts de canvi:
punts_canvi_serie3_ef <-selgmented(serie3_ef.lm, Kmax = 10, type="bic")## No. of breakpoints: 2 .. 3 .. 4 .. 5 .. 6 .. 7 .. 8 .. 9 .. 10 ..
##
## BIC to detect no. of breakpoints:
## 0 1 2 3 4 5 6 7
## 698.6812 701.8040 707.7414 714.1851 721.0583 728.8687 713.0397 710.0976
## 8 9 10
## 714.0422 687.3246 694.5366
##
## No. of selected breakpoints: 8 (1 breakpoint(s) removed due to small slope diff)serie3_ef.seg <- segmented(serie3_ef.lm, seg.Z=~x, psi = c(5,9,16,21,28,35,40,45,54,59))
summary(serie3_ef.seg)$psi## Initial Est. St.Err
## psi1.x 5 5.499227 0.5320942
## psi2.x 9 10.346005 0.5179006
## psi3.x 16 16.491929 0.5806260
## psi4.x 21 23.272538 0.5013678
## psi5.x 28 28.484784 0.4873200
## psi6.x 35 34.969280 0.4593709
## psi7.x 40 40.649312 0.3853447
## psi8.x 45 46.190155 0.3750704
## psi9.x 54 55.999956 0.2048643
## psi10.x 59 58.047110 0.2636145Ens queden els punts de canvi a: 3-2016, 7-2016, 1-2017, 8-2017, 1-2018, 8-2018, 2-2019, 7-2019, 5-2020, 7-2020.
Ara, el valor de \(R^2\) de la segmentada és
summary(serie3_ef.seg)$adj.r.squared## [1] 0.8531528Anem a visualitzar la regressió segmentada damunt les nostres dades
p_serie3_ef <- ggplot(ef[1:63,], aes(x=x, y=y)) + geom_line()+
ylim(c(0,1300))+
ggtitle("Regressió lineal i segmentada sobre les dades \nd'Eivissa i Formentera després de la COVID") +
xlab('índex del mes')+
ylab('Despeses mensuals en €')
my.coef3_ef <- coef(serie3_ef.lm) #coeficients de la recta de regressió lineal
p_serie3_ef <- p_serie3_ef + geom_abline(intercept=my.coef3_ef[1], slope=my.coef3_ef[2], color="green") #afegim la recta de regressió lineal
my.model3_ef <- data.frame(x=1:63, y=fitted(serie3_ef.seg)) #model nou amb els valors ajustats de la regressió segmentada
p_serie3_ef <- p_serie3_ef + geom_line(data=my.model3_ef, aes(x=x,y=y), colour="red") #afegim la recta de regressió segmentada
my.lines3_ef <- serie3_ef.seg$psi[,2]# guardam on estan els punts de canvi estimats
p_serie3_ef <- p_serie3_ef + geom_vline(xintercept=my.lines3_ef, linetype="dashed")
p_serie3_ef <- p_serie3_ef + geom_vline(xintercept = 0)+ geom_hline(yintercept = 0) + theme(panel.background = element_blank()) #eliminam el fons i la quadrícula del gràfic
ggplotly(p_serie3_ef)Per escriure la funció a trossos tenim:
#punts d'intersecció
intercept(serie3_ef.seg)## $x
## Est.
## intercept1 933.81
## intercept2 -308.16
## intercept3 1969.60
## intercept4 -1013.80
## intercept5 4249.90
## intercept6 -2455.40
## intercept7 5932.00
## intercept8 -5595.60
## intercept9 6360.20
## intercept10 -23475.00
## intercept11 5317.20#pendents
slope(serie3_ef.seg)## $x
## Est. St.Err. t value CI(95%).l CI(95%).u
## slope1 -92.589 29.584 -3.1297 -152.340 -32.843
## slope2 133.260 29.584 4.5043 73.509 193.000
## slope3 -86.907 22.363 -3.8861 -132.070 -41.743
## slope4 93.994 17.680 5.3165 58.289 129.700
## slope5 -132.180 29.584 -4.4680 -191.930 -72.436
## slope6 103.220 22.363 4.6155 58.055 148.380
## slope7 -136.630 22.363 -6.1097 -181.800 -91.469
## slope8 146.950 22.363 6.5712 101.790 192.120
## slope9 -111.880 12.078 -9.2638 -136.280 -87.494
## slope10 420.890 66.152 6.3625 287.290 554.480
## slope11 -75.126 29.584 -2.5394 -134.870 -15.380L’error de la regressió segmentada és:
accuracy(serie3_ef.seg)## ME RMSE MAE MPE MAPE MASE
## Training set 9.024968e-16 75.4706 55.09207 -Inf Inf 0.2838393Anem a fer la predicció per a l’any 2021. Recordem que els punts de canvi es donen a 3-2016, 7-2016, 1-2017, 8-2017, 1-2018, 8-2018, 2-2019, 7-2019, 5-2020, 7-2020. Degut a la pertorbació del COVID, sí que hi segueix havent un punt de canvi en estiu i l’altre en hivern successivament però ara no és donen al mateix mes. Per això, per predir els següents punts de canvi agafarem el mes més freqüent. Aleshores els següents punts de canvi seran en 1-2021, 7-2021 i 1-2022.
- La pendent de 7-2020 a 1-2021 i de 7-2021 a 1-2022 serà : -116.899
- La pendent de 1-2021 a 7-2021 serà : 179.663
Els nous punts d’intersecció es calculen de forma anàloga que a Mallorca. Per a les tres noves rectes obtenim que són \(n_1\)=7743.65, \(n_2\)=−11236.32 i \(n_3\)=9523.02.
Aleshores la predicció de 2021 serà:
\(\hat{y}= \left\{ \begin{array}{lcc} -116.899x + 7743.65 & si & x \leq \psi_{11} \\ \\ 179.663x - 11236.32 & si & \psi_{11} < x \leq \psi_{12} \\ \\ -116.899x + 9523.02 & si & \psi_{12} < x \\ \end{array} \right.\)
on \(\psi_7 = 64\) i \(\psi_8 = 70\).
Visualitzem-la:
p3_serie_ef <- ggplot(ef, aes(x=x, y=y)) + geom_line()+
ylim(c(0,1300)) +
ggtitle("Predicció d'Eivissa i Formentera amb el model de \nregressió segmentada després de la COVID") +
xlab('índex del mes')+
ylab('Despeses mensuals en €')
my.model3_ef <- data.frame(x=1:63, y=fitted(serie3_ef.seg)) #model nou amb els valors ajustats de la regressió segmentada
p3_serie_ef <- p3_serie_ef + geom_line(data=my.model3_ef, aes(x=x,y=y), colour="red") #afegim la recta de regressió segmentada
my.lines3_ef <- serie3_ef.seg$psi[,2]# guardam on estan els punts de canvi estimats
p3_serie_ef <- p3_serie_ef + geom_vline(xintercept=my.lines3_ef, linetype="dashed")
p3_serie_ef <- p3_serie_ef + geom_vline(xintercept = 0)+ geom_hline(yintercept = 0) + theme(panel.background = element_blank()) #eliminam el fons i la quadrícula del gràfic
p3_serie_ef <- p3_serie_ef + geom_abline(intercept = 7743.65, slope=-116.899, colour="green") +
geom_abline(intercept = -11236.32 , slope=179.663, colour="blue") +
geom_abline(intercept = 9523.02 , slope=-116.899, colour="orange")
ggplotly(p3_serie_ef)Vegem quin és aquest error que es comet:
o3_ef<-c(serie_ef[64:75]) # dades reals per fer predicció del 2021
v5_ef=c(64:70)
f5_ef <- sapply(v5_ef, function(x) 179.663*x-11236.32) #predicció de gener 2021 a juliol 2021
v6_ef=c(71:75)
f6_ef <- sapply(v6_ef, function(x) -116.899*x + 9523.02) #predicció d'agost 2021 a desembre 2021
p3_ef <- c(f5_ef,f6_ef) #predicció de 2021
sqrt(sum((p3_ef-o3_ef)^2)/13)## [1] 192.499- DECOMPOSE
Visualitzem la sèrie descomposada:
decompose_s3_ef <- decompose(serie3_ef)
plot(decompose_s3_ef, xlab="Any")
Les components de tendència són:
t3_ef <- decompose_s3_ef$trend #tendències de la sèrie sense el tros a predir abans de la COVID
t3_ef## Jan Feb Mar Apr May Jun Jul Aug
## 2015
## 2016 NA NA NA 743.1446 747.8983 758.9400 759.6646 763.9542
## 2017 806.7158 817.0208 824.6025 837.5362 845.8621 843.2208 842.9546 837.8417
## 2018 831.0383 830.7542 836.4517 830.5658 816.6008 804.8000 801.1833 800.7592
## 2019 793.0312 797.2383 795.3000 794.5121 794.6762 798.6983 802.9879 803.0183
## 2020 639.9054 615.4329 594.0237 577.0025 570.0400 569.5333 NA NA
## Sep Oct Nov Dec
## 2015 NA NA NA
## 2016 773.8896 782.0333 791.3829 797.4979
## 2017 833.9504 838.2292 840.2021 837.6517
## 2018 794.4083 787.4583 785.2383 787.7225
## 2019 805.1829 776.5217 710.6350 662.4762
## 2020 NA NA NA NAI les d’estacionalitat:
s3_ef <- decompose_s3_ef$seasonal #estacionalitat
s3_ef <- s3_ef[4:15] #estacionalitat de gener a desembre
s3_ef## [1] -226.71583 -201.30468 -200.48260 -185.77737 -83.92262 118.54038
## [7] 314.53677 311.05855 156.86407 185.77125 -73.42270 -115.14520Anem a fer la predicció:
a3_ef <- c(s3_ef[7:12],s3_ef) #estacionalitat de juliol 2020 a desembre 2021 (es per poder fer la predicció)
pred3_decompose_ef <- sapply(a3_ef, function(x) 569.5333 + x) #predicció de juliol de 2019 a desembre 2020
p3_dec_ef <-c(serie3_ef[1:57], pred3_decompose_ef)
A3_ef<- data.frame("x" = ef[1:75,]$x, "y" = ef[1:75,]$y, "p"= p3_dec_ef)
grafica3_ef_dec <- ggplot(A3_ef)+
geom_line(aes(x,p), color="red")+
geom_line(aes(x,y))+
labs(title="Predicció d'Eivissa i Formentera després de la COVID amb el model \nde descomposició", x="Índex del mes", y="Despeses mensualns en €")+
geom_vline(xintercept=0) + geom_hline(yintercept=0) + theme(panel.background = element_blank())
grafica3_ef_dec
Calculem lerror de la predicció
sqrt(sum((c(serie_ef[64:75]- pred3_decompose_ef[7:18]))^2)/12)## [1] 256.8545Amb la funció predict() la predicció seria la
següent:
prediccio3_ef <- predict(serie3_ef,h=12)
plot(prediccio3_ef, xlab="Any", ylab="Despeses mensuals en €")
Vegem-la juntament amb les nostres dades:
df_prediccio3_ef <- data.frame(prediccio3_ef)
p3_ets_ef <- data.frame("x"= 1:75, "PointForecast"=c(serie3_ef,df_prediccio3_ef$Point.Forecast), "Lo80" = c(rep(NA,63),df_prediccio3_ef$Lo.80), "Hi80" = c(rep(NA,63),df_prediccio3_ef$Hi.80), "Lo95" = c(rep(NA,63),df_prediccio3_ef$Lo.95),"Hi95" = c(rep(NA,63),df_prediccio3_ef$Hi.95))
grafica_pred3_ets <- ggplot((ef[1:75,]))+
geom_ribbon(data = p3_ets_ef, aes(x, ymin = Lo95, ymax = Hi95), fill = "blue", alpha = 0.25) +
geom_ribbon(data =p3_ets_ef, aes(x, ymin = Lo80, ymax = Hi80), fill = "blue", alpha = 0.25) +
geom_line(data = p3_ets_ef, aes(x, PointForecast), colour = "blue") +
geom_line(aes(x,y), color="red")+
labs(title="Predicció del model ETS(A,N,A) a Eivissa i Formentera \ndesprés de la COVID", x="Índex del mes", y="Despeses mensuals en €")+
geom_vline(xintercept=0) + geom_hline(yintercept=0) + theme(panel.background = element_blank())
grafica_pred3_ets
Obtenim un error d’uns 138 euros.
accuracy(prediccio3_ef,serie_ef)## ME RMSE MAE MPE MAPE MASE ACF1
## Training set -2.165247 122.8461 85.26169 -Inf Inf 0.6811543 0.1330175
## Test set 115.747684 137.8193 123.36967 15.6055 16.37964 0.9855984 0.5931860
## Theil's U
## Training set NA
## Test set 1.944982- SARIMA
Quin model proposam nosaltres?Vegem l’ACF i el PACF:
par(mfrow=c(1,2))
acf(serie3_ef)
pacf(serie3_ef)
Per a la part regular: p=2, q=2 i d=0.
Sí hi ha indicació d’estacionalitat, llavors feim una diferenciació d’ordre 12 i tornam a calcular l’ACF i el PACF:
serie3_ef_diff <- diff(serie3_ef,12)
plot.ts(serie3_ef_diff, main="Sèrie sense estacionalitat", ylab="Sèrie diferenciada", xlab="Any")
par(mfrow=c(1,2))
acf(serie3_ef_diff)
pacf(serie3_ef_diff)
Aleshores obtenim un model ARIMA(2,0,2)(1,1,3)[12]
R proposa el següent model:
auto.arima(serie3_ef)## Series: serie3_ef
## ARIMA(0,1,2)(0,1,1)[12]
##
## Coefficients:
## ma1 ma2 sma1
## -0.1957 -0.5346 -0.7096
## s.e. 0.1395 0.1384 0.3900
##
## sigma^2 = 20382: log likelihood = -322
## AIC=651.99 AICc=652.88 BIC=659.64Llavors tenim els següents models
model5_ef <- arima(serie3_ef, order=c(2,0,2), seasonal = list( order=c(1,1,3), period=12))
model6_ef <- arima(serie3_ef, order=c(0,1,2), seasonal = list( order=c(0,1,1), period=12))Amb uns errors de:
accuracy(model5_ef)## ME RMSE MAE MPE MAPE MASE ACF1
## Training set -7.271777 109.7292 64.51072 -Inf Inf 0.4889004 0.01064585accuracy(model6_ef)## ME RMSE MAE MPE MAPE MASE ACF1
## Training set -13.9712 123.312 71.6025 -Inf Inf 0.5426461 -0.03379418I les seves prediccions són:
fc_m5_ef <- forecast(model5_ef, h=12)
fc_m6_ef <- forecast(model6_ef,h=12)#primer model
pre3_ef <- data.frame("Point Forecast" = serie3_ef, "Lo 80" = rep(NA,63), "Hi 80"= rep(NA,63), "Lo 95" = rep(NA,63), "Hi 95" = rep(NA,63))
pred_m5_ef <-data.frame(fc_m5_ef)
grafica_m5_ef <- data.frame("x" = 1:75, "PointForecast" = c(pre3_ef$Point.Forecast,pred_m5_ef$Point.Forecast), "Lo80" = c(pre3_ef$Lo.80, pred_m5_ef$Lo.80), "Hi80"= c(pre3_ef$Hi.80, pred_m5_ef$Hi.80), "Lo95" = c(pre3_ef$Lo.95, pred_m5_ef$Lo.95), "Hi95" = c(pre3_ef$Hi.95, pred_m5_ef$Hi.95))
grafica5_ef <- ggplot(ef[1:75,])+
geom_ribbon(data = grafica_m5_ef, aes(x, ymin = Lo95, ymax = Hi95), fill = "blue", alpha = 0.25) +
geom_ribbon(data = grafica_m5_ef, aes(x, ymin = Lo80, ymax = Hi80), fill = "blue", alpha = 0.25) +
geom_line(data = grafica_m5_ef, aes(x, PointForecast), colour = "blue") +
geom_line(aes(x,y), color="red")+
labs(title="Predicció d'Eivissa i Formentera amb el model ARIMA (2,0,2)(1,1,3)[12] \ndesprés de la COVID", x="Índex del mes", y="Despeses mensuals en €")+
geom_vline(xintercept=0) + geom_hline(yintercept=0) + theme(panel.background = element_blank())
grafica5_ef
#segon model
pred_m6_ef <-data.frame(fc_m6_ef)
grafica_m6_ef <- data.frame("x" = 1:75, "PointForecast" = c(pre3_ef$Point.Forecast,pred_m6_ef$Point.Forecast), "Lo80" = c(pre3_ef$Lo.80, pred_m6_ef$Lo.80), "Hi80"= c(pre3_ef$Hi.80, pred_m6_ef$Hi.80), "Lo95" = c(pre3_ef$Lo.95, pred_m6_ef$Lo.95), "Hi95" = c(pre3_ef$Hi.95, pred_m6_ef$Hi.95))
grafica6_ef <- ggplot(ef[1:75,])+
geom_ribbon(data = grafica_m6_ef, aes(x, ymin = Lo95, ymax = Hi95), fill = "blue", alpha = 0.25) +
geom_ribbon(data = grafica_m6_ef, aes(x, ymin = Lo80, ymax = Hi80), fill = "blue", alpha = 0.25) +
geom_line(data = grafica_m6_ef, aes(x, PointForecast), colour = "blue") +
geom_line(aes(x,y), color="red")+
labs(title="Predicció d'Eivissa i Formentera amb el model ARIMA (0,1,2)(0,1,1)[12] \ndesprés de la COVID", x="Índex del mes", y="Despeses mensuals en €")+
geom_vline(xintercept=0) + geom_hline(yintercept=0) + theme(panel.background = element_blank())
grafica6_ef
Vegem quin és l’error de cada model i estudiem-los:
accuracy(fc_m5_ef,serie_ef[64:75], h=12)## ME RMSE MAE MPE MAPE MASE
## Training set -7.271777 109.7292 64.51072 -Inf Inf 0.4889004
## Test set 129.020379 164.9598 136.90003 17.47361 18.27666 1.0375094
## ACF1
## Training set 0.01064585
## Test set NAaccuracy(fc_m6_ef,serie_ef[64:75], h=12)## ME RMSE MAE MPE MAPE MASE ACF1
## Training set -13.9712 123.3120 71.6025 -Inf Inf 0.5426461 -0.03379418
## Test set 255.6337 270.3203 255.6337 33.14161 33.14161 1.9373431 NAcheckresiduals(fc_m5_ef)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(2,0,2)(1,1,3)[12]
## Q* = 3.0367, df = 5, p-value = 0.6943
##
## Model df: 8. Total lags used: 13checkresiduals(fc_m6_ef)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(0,1,2)(0,1,1)[12]
## Q* = 5.8462, df = 10, p-value = 0.828
##
## Model df: 3. Total lags used: 13par(mfrow=c(1,2))
qqPlot(fc_m5_ef$residuals, id=FALSE, mean=mean(fc_m5_ef$residuals), sd=sd(fc_m5_ef$residuals))
qqPlot(fc_m6_ef$residuals, id=FALSE, mean=mean(fc_m6_ef$residuals), sd=sd(fc_m6_ef$residuals))
shapiro.test(fc_m5_ef$residuals)##
## Shapiro-Wilk normality test
##
## data: fc_m5_ef$residuals
## W = 0.77159, p-value = 1.605e-08shapiro.test(fc_m6_ef$residuals)##
## Shapiro-Wilk normality test
##
## data: fc_m6_ef$residuals
## W = 0.78424, p-value = 3.217e-08lillie.test(fc_m5_ef$residuals)##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: fc_m5_ef$residuals
## D = 0.15381, p-value = 0.0007865lillie.test(fc_m6_ef$residuals)##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: fc_m6_ef$residuals
## D = 0.16203, p-value = 0.0002898Resum dels errors que comet cada un del models anteriors:
| reg. segmentada | descomposició clàssica | ETS(A,N,A) | ARIMA (2,0,2)(1,1,3)[12] | ARIMA (0,1,2)(0,1,1)[12] | |
|---|---|---|---|---|---|
| error model | 75.47 | 122.85 | 109.73 | 123.31 | |
| error predicció | 192.499 | 348.405 | 137.82 | 164.96 | 270.32 |